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Abstract 

We investigate the thermalization of charm quarks in high energy heavy ion colhsions. To this 
end, we calculate the diffusion coefficient in the perturbative Quark Gluon Plasma and relate it 
to collisional energy loss and momentum broadening. We then use these transport properties to 
formulate a Langevin model for the evolution of the heavy quark spectrum in the hot medium. 

— 1/2 

The model is strictly valid in the non-relativistic limit and for all velocities jv < as to leading 
logarithm in T/m^- The corresponding Fokker-Planck equation can be solved analytically for a 
Bjorken expansion and the solution gives a simple estimate for the medium modifications of the 
heavy quark spectrum as a function of the diffusion coefficient. Finally we solve the Langevin 
equations numerically in a hydrodynamic simulation of the heavy ion reaction. The results of this 
simulation are the medium modifications of the charm spectrum Raa and the expected elliptic 
flow V2{pt) as a function of the diffusion coefficient. 



1 



I. INTRODUCTION 



The experimental relativistic heavy ion program at RHIC aims to measure the properties 
of the Quark Gluon Plasma (QGP) [1]. One of the most exciting results from RHIC so far is 
the large azimuthal anisotropy of light hadrons with respect to the reaction plane, known as 
elliptic flow. Elliptic flow has been measured as a function of impact parameter, transverse 
momentum, rapidity and particle type and is quantified with V2{pt) 



/ d^j£^, cos(20) 



V2{Pt) = ' ^^'^''^Zn^^'^' ■ (1-1) 



The observed elliptic flow is significantly larger than was originally expected from kinetic 
calculations of quarks and gluons L Zl , but in fairly good agreement with simulations based 
upon ideal hydrodynamics JSI, 0, llClL . This result suggests that the medium responds as 
a thermalized fluid and that the transport mean free path is small 0, 0, 0| . 

However, this interpretation of the elliptic flow results is not universally accepted 0, 
isl - Hadronization may amplify the underlying partonic elliptic flow Indeed, parton 



coalescence is one mechanism which may amplify the hadronic elliptic flow relative to its 
partonic constituents [H, 17, [3, [l^ 2^. On physical grounds, it seems unlikely that the 



typical mean free path is much smaller than a thermal wavelength, 1 / (27rT) . Indeed it has 
recently been conjectured that the hydrodynamic diffusion coefficient rj/{e + p) is strictly 



larger than half a thermal wavelength j21 



e + p 47rT 

Here r] is shear viscosity, (e + p) is the enthalpy, and the ratio is a fundamental length scale 
in the QGP. \i ri/{e + p) is significantly larger than this conjectured bound hydrodynamics 
would not be a viable explanation for the observed flow at RHIC. In this work we will accept 
the hydrodynamic paradigm of the RHIC results and study the correlated consequences of 
this interpretation. 

Heavy quarks are a good probe of the transport properties of the medium. Given an 
estimate of the light quark relaxation time ^ r]/{e + p), the heavy quark relaxation time is 

M 7] 

tr 



T e + p 



where M is the mass of the heavy quark and T is the temperature. Thus, with M ^ 1.4 GeV 
and T ^ 250 MeV, we expect that the charm equilibration time is approximately 6 times 
larger than the light quark equilibration time. Since this factor is relatively large, we further 
expect that the elliptic flow of charm quarks will be smaller than the flow of light hadrons. 

In addition, the heavy quarks are produced with a power-law transverse momentum 
spectrum which deviates strongly from the thermal spectrum. The relaxation time will 
control the extent to which the initial power-law spectrum approaches the thermal spectrum. 
Similarly, the relaxation time will control the extent to which the charm quark will follow 
the underlying flow of the medium. If the charm quark completely follows the flow of the 
medium then thermal spectrum is actually quite close to the perturbative spectrum and 
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it may be difficult to distinguish these two cases Medium modifications of the 

heavy quark spectrum Raa flow will be studied ex p eriment ally this year and will provide 
an experimental estimate of this relaxation time [23|, |24 , Hill 113. 

Most recent studies of the medium modiflcations of the charm spectrum have computed 
the energy loss of a heavy quark by gluon bremsstrahlung |31|] . In weak coupling 

(which is the framework in which all calculations have been performed), bremsstrahlung is 
the dominant energy loss mechanism only if the heavy quark is ultra-relativistic, 7^ ^ 1/ g. 
(Similarly, for an electron traversing a hydrogen target, bremsstrahlung losses flrst exceed 
ionization losses when 7^ ~ 700 I32l].) For much of the measured momentum range, the 
heavy quark is not ultra-relativistic, 7^ < 4, and in this case it is far from clear that radiative 
energy loss dominates over collisional energy loss. 

About two thirds of all heavy quarks are produced with p < M, and therefore radiative 
energy loss should be neglected when studying bulk thermalization. When 'yy > 4, calcula- 



tions do suggest that radiation dominates the average energy loss rate [33,|3J|. However, as 
has been repeatedly emphasized (sl,!!^, the average energy loss is insufficient to describe 
the medium modiflcations of the spectrum Raa- Collisions have a different fluctuation spec- 
trum than radiation and therefore might contribute more to the suppression factor than was 
at first anticipated Since we are primarily interested here in heavy quarks with typical 
momenta 71; ~ 1, we will concentrate exclusively on elastic collisions. 

Considering these points, we will re-examine collisional energy loss of a heavy quark in 
the perturbative QGP. Our tools are Hard Thermal Loop (HTL) perturbation theory and 
a heavy quark expansion (M ^ T). The average energy loss rate was first computed by 
Braaten and Thoma [s^ and we will independently verify their results. (Recently this cal- 
culation was extended to anisotropic plasmas by Romatschke and Strickland 0].) We will 
also compute the rates of longitudinal and transverse momentum broadening which are es- 
sential to a complete calculation of the modification factor Raa- We will relate all of these 
rates to the diffusion coefficient which we will compute. (In principle the diffusion coefficient 



of a heavy quark could have been gleaned from the results of Braaten and Thoma j38| and 
Svetitsky j40].) In perturbation theory, we can compare the diffusion coefficient of the heavy 
quark to the hydrodynamic time scale ?7/(e + p) which was calculated previously 41, 4^ 43| . 
Many of the ambiguities of perturbation theory cancel in the ratio of transport coefficients 
and we therefore hope to be able to extrapolate smoothly into the non-perturbative do- 
main. Following this ideology, we express all of our phenomenological results in terms of the 
diffusion coefficient which may ultimately be determined from lattice QCD calculations. 

With these transport properties in hand, we adopt a Langevin model for the equilibration 
of heavy quarks in heavy ion collisions. The Langevin equations correctly describe the 
kinetics of a heavy particle in a thermal medium and therefore naturally interpolate between 
a hydrodynamic regime at small momentum and a kinetic regime at large momentum. The 
model is similar to old work by Svetitsky ji^ and was later used without fluctuations 
to estimate Raa Js^]. The model is strictly valid for non-relativistic quarks and for all 
velocities to leading logarithm mT /mo- The corresponding Fokker-Planck equation is solved 
analytically in Section |3 for a Bjorken expansion. The solution provides a simple estimate 
for the modification factor Raa and further elucidates the dynamics of equilibration. Then 
we solve the Langevin equations numerically in a hydrodynamic simulation of the heavy ion 
reaction. The results of the simulation are the medium modification factor Raa and the 
corresponding elliptic fiow 1^2 (pt) as a function of the diffusion coefficient. 
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Very recen tly, two papers have adopted a similar Langevin approach to address charm 
equilibration ji^, H^- The first of these papers 44 1 estimated possible non-perturbative 
contributions to the drag coefficient which may arise from quasi-hadronic bound states 
in the QGP. The second of these papers estimated radiative and collisional energy 
loss, and found that collisional loss is significant even for rather energetic charm quarks, 
E ~ 5 - 10 GeV. 

In addition, classical Boltzmann simulations by Molnar ji^ and subsequent simulations 
by Zhang et al. ji^] have studied how V2{pt) and studied how V2{pt) and Raa depend on 
the charm mean free path. As is discussed in Section FV C[ the results of Molnar's simulation 
are comparable to the Boltzmann-Langevin approach adopted here. 

Throughout, we will denote 4- vectors with capital letters P, Q and use p, q for their 
3- vector components, p^,q^ for their energy components, and p,q for |p|, |q|. Our metric 
convention is [-,+,+,+]. 



II. NON-RELATIVISTIC HEAVY QUARKS IN A THERMAL MEDIUM 

First consider thermal heavy quarks, M ^ T, with typical thermal momentum p ~ 
\fMT and velocity v ~ \JT/M <C 1. Since p ^ T it takes many collisions to change 
the momentum substantially. Even for hard collisions with momentum transfer g ~ T, it 
takes ~ M/T collisions to change the momentum by a factor of order one. Therefore it 
is a good approximation to model the interaction of the heavy quark with the medium as 
uncorrelated momentum kicks. The momentum of the heavy quark will evolve according to 
the macroscopic Langevin equations ji^ 



^ = m - VDPr , mMt')) = K5,,5{t - t') . (2.1) 

Here rjo is a momentum drag coefficient and ^i{t) delivers random momentum kicks which 
are uncorrelated in time. 3k is the mean squared momentum transfer per unit time. The 
solution of this stochastic differential equation is 



p,{t) = / dt'e^-^''~'^Ut') , (2.2) 

J —oo 

where we have assumed that t ^ i]^^. The mean squared value of p is 

3MT = ip') = f dhdhe^-^'^+'-^\Uti)Ut2)) = P- , (2.3) 

and therefore 

Vd = — ^ • (2.4) 

Now the diffusion constant in space, D, can be found by starting a particle at a; = at 
t = and finding the mean squared position at a later time, 

{xi{t)xj{t)) = 2DtSij QDt = (x^(t)) . (2.5) 
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Using the relation between position and momentum 

M 



x,{t)= I'dt'P^, (2.6) 







we have 

ft rt 



6a = ^*^<ife-jlj{p(t,)p(*,)) = -^. (2.7) 

and therefore the diffusion constant is 0| 

T 2T^ 

D = T^ = — - (2.8) 

In the next paragraphs we will determine the mean squared momentum transfer per unit 
time, Sk, and then the diffusion coefficient. 

The only 2 2 scattering processes are qH qH {H the heavy quark) and gH gH. 
The qH qH process only occurs by t channel gluon exchange. Since in the rest frame of the 
plasma the Compton amplitude is suppressed by w ^ ~ T/M (see Appendix|B)), the gH gH 
process is also dominated by t channel gluon exchange. Kinematics demand that the transfer 
momentum be spatial up to 0{v) corrections, so the Hard Thermal Loop correction on the 
transferred gluon is described simply by Debye screening. Writing the incoming and outgoing 
momenta of the thermal particle as k and k', and taking that particle's dispersion relation 
to be ultra-relativistic, the scattering matrix elements squared, summed on the colors and 
spins of the incoming thermal particles and over all quantum numbers of the final state 
particles, are (see Appendix rB|). 



\M\ 



2 

quark 



Ch9 



IGM^kl (1 + cos ^kk' 



(g^ + ^'^ 



D) 



l-^lgiuon = [NcCng'] IQM'kl (1 + cos^ gickO /^2^)2 ' (2-9) 

Ch = C*F denotes the color Casimir of the heavy quark and an extra factor of 2 for the quark 
case accounts for anti-quarks. 

The mean squared momentum transfer per unit time is 3k. To compute this quantity the 
matrix elements must be integrated over the incoming momentum k and outgoing momenta 
k', p', weighted by the appropriate statistical functions and by the squared momentum 
transfer = (k — k')^. The mean squared momentum transfer per unit time (3k) is then 

1 f rl^\<:rl^\<:'rl^W 

'^-^I il^m? (^-'^^^'p - - p' - "'^'^f*' - ^ 

X [A'/|A^I?u„kn/(*)(l-n/(«:')) + |A1|J„„r!t(A)(l+n4(<:'))] . (2.10) 

In writing this formula we have used the non-relativistic limit = = M. It is convenient 
to shift the p' integration to an integral over the momentum transfer q = p' — p ; the 
momentum conserving delta function becomes 6^(k' — q — k). The appendix shows how a 
simple change of variables and an expansion in <^ makes these integrals relatively 
straightforward. Using the relation between k, and the diffusion coefficient Eq. ()2.8j) . we find 



D 



ChQ'T 



, 2T 1 ^C'(2)\ ATf/ AT 1 C{2)\ 

In + 77 - 7e + -TTTy- +— In + tt-7e + 777^ 

mo 2 C(2) / 2 V "^D 2 C(2) / 



f2.111 
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FIG. 1: (Color online) (a) The diffusion coefficient of a heavy quark in a QGP with Nf flavors 
of light quarks, (b) The ratio of the diffusion coefficient of a heavy quark to the hydrodynamic 
diffusion coefficient r]/{e + p). 



This formula for the diffusion coefficient could have been extracted by combining the cal- 
culations of Braaten and Thoma js^ and Svetitsky ji^ . Corrections to this expression are 
suppressed by at least one power of y/ct^; we expect corrections at that level. 

This expression for the diffusion coefficient is based on a small mo expansion and therefore 
becomes unreliable when the Debye mass becomes large. Rather than using a small m-D 
expansion we can numerically integrate Eq. ()2.10j) to determine the the diffusion coefficient. 
Something of this sort is necessary to deal with large values of mD, but we emphasize that, 
at large values of m^, the procedure is ad-hoc and should be considered only a quahtative 
guide. The mD expansion would be a good approximation if the coupling were truly small. 
The numerical evaluation of the diffusion coefficient is illustrated in Fig. ^a) as a function 
of mD. 

As discussed in the introduction, the time scale for equilibration r ~ Y'^' ^o'^ l^t us 
make this more concrete. The rate of equilibration is given by 77^^, as can be intuited from 
the Langevin equations, Eq. ()2.1|1 . or from the analysis of Sectional Vd^ is directly related 
to the diffusion coefficient via 

1 M , , 

— = ^D. 2.12 

Taking as ~ 0.5, rriD/T ^ 1.5, and M/T ^ 7 we estimate that the diffusion coefficient is 
D ^ ^ ^ The relaxation time is then r^^^ ^ Q.7 /T. 

It is useful to compare this timescale with other hydrodynamic relaxation times in the 
QGP. The damping rate of sound waves in the QGP is controlled by f]/{e + p) where rj is 
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the shear viscosity and e + p is the enthalpy. Dividing by this ratio we have 



1 



M 7] 



D 



(2.13) 



X 



r]/{e + p) 



T e + p 



The quantity in square brackets is an estimate of the ratio between the heavy quark and 
hydrodynamic relaxation times and is illustrated in Fig. ^b) as a function of mo/T. Here 
the shear viscosity at leading order has been taken from 43]. As seen from the figure, for 
Nf ^ 3 the diffusion coefficient is ~ 6 times larger than hydrodynamic scale ri/{e + p). 

Both the calculation of the diffusion coefficient D and of the shear viscosity t] are plagued 
by ambiguities when the coupling is strong, as > 0.2. In particular, in both processes it 
is unclear how to screen the plasma when ttid is large and how to estimate the form and 
size of subleading corrections. These issues have been discussed in |4j|, but were not 
resolved. However, since both processes are dominated by t-channel gluon exchange, the 
ambiguities in the calculation should largely cancel in the ratio of transport coefficients. We 
therefore hope that the computed factor of ~ 6 for D/[r]/{e+p)] is largely independent of its 
perturbative assumptions. Thus, with ^ ^ 7, T ^ 200 MeV and an optimistic estimate of 
the shear viscosity, ^ = 1/(6T), we estimate that the charm quark thermalization time is 
of order ^ 7 fm. This time scale should be compared with the time scale for the development 
of elliptic flow ^ 4 fm. 

III. ENERGY LOSS 

Since the initial distribution of charm quarks is much "harder" than a thermalized spec- 
trum, it is important to study how higher energy heavy quarks, with 7^ ~ 1, lose their 
energy in the thermal medium. We now turn to a discussion of this problem. Our discussion 
differs from most recent literature in that we take the dominant energy loss mechanism to 
be elastic scattering. 

A. Why 2 <-> 2 dominates for 71; ~ 1 

When the coupling is small, bremsstrahlung dominates the energy loss rate for very 
fast particles (with jv ^ 1/g) while collisions dominate the rate for moderately relativistic 
particles with 7?; ~ 1. 

To see why, consider a heavy particle of momentum v, undergoing a scattering. Call its 
momentum P, with p = |p| = vp^. For energy loss in a condensed medium, the scattering 
is typically off a nucleus, so the transferred 4-momentum Q^^ is purely spatial. In a thermal 
medium like the QGP, the scattering is off a relativistic particle at a random angle with 
respect to the heavy particle; the kinematics of the other particle requires that be 
spacelike, but can be ~ g = |q|. Kinematics requires that the lightlike, bremmed particle's 
momentum K satisfy {K + P') = {Q + P), with P' the final particle momentum. K can be 
the largest if it is perfectly coUinear with P and if Q is anti-coUinear; in this case, energy 
conservation reads = p^ — p'^ + = v{p — p') + , while momentum conservation reads 
k = = p — p' — q. For the spatial (Coulomb) scattering case, k < vq/{l—v), while for the 
QGP case, k < {vq+q^)/{l—v), so the actual energy loss is p° — < v{q+q^)/{l—v). In 
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either case, for 7^ ~ 1, the energy of a bremmed particle is at best of order of the transfer 
momentum in the scattering. 

The energy loss in a scattering process without bremsstrahlung is When g° ~ T the 
plasma temperature, phase space favors events with negative and of order q. Since the 
rate of bremsstrahlung is suppressed by an additional factor of ctg, and the typical energy 
loss is the same, the contribution of bremsstrahlung to energy loss is smaller by a factor of 
ctg. In the soft region, q ^ gT, there is an approximate cancellation between energy losing 
and energy gaining scatterings, and the dominance of 2 ^ 2 processes is only by v^g, but 
they remain dominant. Softer momentum transfers are screened by the plasma and do not 
play an important role in energy loss. 

Very similar arguments apply for electromagnetic energy loss of a high energy particle in 
a medium, except that the scale gT is played by the atomic radius and the 2 2 processes 
do not have any cancellation between energy losing and gaining processes; for 7?; ~ 1, 
bremsstrahlung energy losses are subdominant to ionization losses (due to elastic scattering 
with electrons, transferring more energy than the electron's binding energy) by Za. 

As the particle's energy increases, the importance of bremsstrahlung also increases. While 
the bremsstrahlung cross-section remains suppressed by a power of ctg, the amount of energy 
which can be lost in a single event increases. As we will see, the rate of energy loss by 2 2 
processes only increases logarithmically in 'jv. Bremsstrahlung energy losses typically rise 
almost linearly in jv, and so become dominant at ■jv ~ 1/Za for electromagnetic processes 
and at 71) ~ 1/ ^/a^ for heavy quark energy loss in the QGP. Since only a small tail of heavy 
quarks have ^yv 1, we will only treat 2 <-> 2 processes in this work. 



B. Relativistic heavy quarks; momentum loss and momentum diffusion 

First consider a heavy quark in a static medium with p ^ T and velocity 7f ~ 1. It 
takes ~ p/T collisions to change the momentum of the heavy quark by a factor of order one. 
The time between hard collisions is ~ l/{g'^T) and thus the equilibration time scale a heavy 
quark is of order ~ {p/T) l/{g*T). 

Next consider how a heavy quark with momentum p ^ T changes over a time interval 
At which is long compared to the timescale of medium correlations but short compared to 
the time scale of heavy quark thermalization, 

« At « ^— - . (3.1) 



g4rp rp g4rp " 

For At large compared to l/{g'^T) the number of collisions is large ~ {At){g'^T). On 
the other hand the momentum of the heavy particle has scarcely changed since the total 
momentum transfered Ap is of order Ap ~ T[At){g^T) which is small compared to p. This 
means that the probability distribution for a given momentum transfer from a single collision 
is approximately constant over this time period. The accumulated momentum transfer is a 
sum of A^ such collisions. The sum of a large number of momentum transfers drawn from 
an identical probability distribution is approximately Gaussian plus corrections which go 
as 1/A^. In the next At time interval the process repeats itself independently. Thus we 
can write down a macroscopic equation of motion for the the heavy quark moving in the z 
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direction, 



— (p) = -r]D{p)p, 

Here ((Apy)^) = is the variance of the momentum distribution transverse to the direc- 
tion of the heavy quark, ((Ap^)^) = (pz — (p^))^ is the variance of the momentum distribution 
in the direction parallel to the direction of the quark, and the time derivatives are under- 
stood to act only on a time scale of order At. The factor of | in the transverse fluctuations 
has been inserted because there are two perpendicular directions. The functions r^/j, kt, and 
Kl encode the average momentum loss and the transverse and longitudinal fluctuations. 

Now we will compute these coefficients using kinetic theory. First we will compute the 
mean rate of momentum loss. This is most easily done by computing the energy loss rate, 
dp^/dt, which is related to the momentum loss by dp^/dt — vdp/dt. The energy loss rate is 
found by multiplying the scattering rate with transfer energy q^, schematically, 

^ = - / \MWm{l±f[k-q']). (3.2) 

dt V Jk,q 

The factor g° is not enough to render this IR convergent, but there are cancellations be- 
tween > Q and g° < contributions. It is easiest to account for these by averaging in 
the integrand over the process with k incoming and k' outgoing, and the process with k' 
incoming, k outgoing, and opposite Q^. Because we take p ^ k the kinematics and matrix 
element are the same, but the population functions differ, yielding 

tt^Yvj,, '-^''^^ ^^^^^^^ ^ " ± m)] ■ (3.3) 

The population function here is (e^/-^ — e^^~'^°^/^) f[k\f[k — which vanishes at small q^ 
and makes the integral well behaved. 

Similarly, to compute the rate of transverse momentum broadening, we weight the tran- 
sition rate with the square of the transverse momentum transfer 

I <(Ap^)^> = J^jM\Vrf[m ± f[k-q']) ■ (3.4) 

Here the integral will be convergent and we do not need to symmetrize over forward and 
backward collisions. Finally the rate of longitudinal momentum broadening is 

I <(Ap,)2> = j^\MUf[k]{l ± f[k-q']) . (3.5) 

Again this integral is convergent. Thus to compute the transport coefficients r]D^ k,t, i^l we 
need to specify the matrix elements and perform (numerically) the phase space integrals. 
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The covariant expressions for the scattering matrix elements in Eq. ()2.9|) (still summed 
over spins and colors of the bath particle), in vacuum, are 



I I quark 





16 


2 





\M\' 



gluon 



[N^ChQ^] 16 



M' 



2 H 



2g2 



+ 



(3.6) 



where Ch = C*f is the quadratic Casimir of the heavy quark, P is the heavy particle 4- 
momentum, K is the light particle 4-momentum, Q is the 4-momentum transfer, and we 
use [-,+,+,+] metric convention. The inclusion of Hard Thermal Loops is only necessary at 
small Q^, where the leading term dominates; the details appear in the appendix. This 
introduces the Debye mass into the problem and thus the results will generally depend on 
the ratio oimo/T. 

Analytic expressions for the transport coefficients can be derived to leading logarithm 
in T /uir). Without the Hard Thermal Loop correction the phase space integrals are loga- 
rithmically divergent at small q^. The leading log transport coefficient is found by evaluat- 
ing the contribution to the phase-space integrals from the logarithmically divergent region 
-C -C T^. This integral will give a log of T/mj) times some function of v. Further de- 
tails are given in the appendix. We find that, to leading logarithm, the transport coefficients 
are 



dp 
'dt 



Nc + — 
2 



2Tdp 
V dt 



N 



247r 
127r 



2f3 



In 



1+f 



1 

2^ 



-V 



4^3 



1-f 

In - 



ln(r/mD) , 

ln(r/mD) 



(3.7) 
(3.8) 
(3.9) 



We note that the relation = is what is expected of the Boltzmann-Langevin ap- 
proach. Indeed, the conclusion of the next section is that the Boltzmann-Langevin approach 
is strictly valid only in a leading-log approximation or in the non-relativistic limit. For com- 
pleteness, the constant under the logarithm in Eqs. ()3.7|) - ()3.9|) is calculated in the appendix. 

As for the diffusion coefficient, these expressions are only valid when the Debye mass 
is very small compared to the temperature. To extrapolate to finite m-o/T we return to 
the original expressions and numerically perform the phase-space integrals to determine 
the transport coefficients. The phase space integrals can be reduced to three dimensional 
integrals as is detailed in the appendix. The resulting transport coefficients and their de- 
pendence on the heavy quark momentum are illustrated in Fig. |21 and are discussed below. 



The momentum dependence of the drag coefficient for various values of mo/T is shown 
in Fig. Efa). As above, we emphasize that these results are strictly valid only when the 
Debye mass is small and therefore the different curves illustrate the uncertainty in the 
calculation. The drag coefficient fjoip) has units (time)~^ and sets the equilibration rate. 
Since the drag coefficient at zero momentum is related to the diffusion coefficient which 
has already been discussed, we have divided by the drag coefficient at zero momentum 
to isolate the momentum dependence. We see that the equilibration rate rjD{p) does not 
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0.5 1 1.5 2 2.5 3 0.5 1 1.5 2 2.5 3 

p/M P/M 

FIG. 2: (Color online) (a) The drag coefficient as a function of the heavy quark momentum, 
^ = floip)?- (b) The transverse {kt{p)) and longitudinal {kl{p)) momentum diffusion coefficients 
as a function of heavy quark momentum for m£,/T = 1.5 . For comparison we also show the 
longitudinal momentum diffusion coefficient in a Boltzmann-Langevin approach 



decrease significantly with momentum. Only when the momentum is larger than 3M does 
the drag coefficient decrease. For comparison we have illustrated the expected form of the 
drag coefficient for an extreme model which is studied later. In this model, we have ^ oc w 
and therefore the equilibration time is inversely proportional to the energy of the heavy 
quark, ?7£) oc 1/E. 

Next we report on the transverse and longitudinal fluctuations for the heavy quark in 
Fig. |21^b). We see that both the transverse and longitudinal fluctuations (kt(p)/^t(0) and 
i^l{p) / rise with the momentum of the heavy quark. It is instructive to compare 
the longitudinal fluctuations derived from the full Boltzmann equation to the longitudinal 
fluctuations expected in a Boltzmann-Langevin approach described in the next section. In 
the Boltzmann-Langevin approach the longitudinal fluctuations are directly related to the 
drag through the relation 

d x2\ 2T dp 

We see that the Boltzmann-Langevin prediction for the longitudinal fluctuations ^) 
underestimates the fluctuations of the full Boltzmann equation (kl)- At large momenta the 
underestimate can be a factor of 2 to 4. 



IV. THE BOLTZMANN FOKKER-PLANCK EQUATION 

The discussion in the previous sections (especially the first paragraph of Section IIII B|l 
suggests a relativistic generalization of the non-relativistic Langevin equations. To this end. 
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we write a stochastic equation of motion for the heavy quark in the rest frame of the medium 

dpL 



dt 
dpr 
dt 



-VdW + ^l, (4.1) 
. (4.2) 



Here dpi and dpT are the momentum increments parallel and transverse to the direction of 
the heavy quark respectively, and are random momentum kicks in the longitudinal 
and transverse directions which satisfy 

(ei(t)ei(O) = KL{p)fiP6{t-t'), (4.3) 

{iW) iW) = ^t{p) {S'^ - mS{t - t') , (4.4) 

mt)eL{t')) = 0. (4.5) 

Thus r]j:>{p)p is the momentum loss per unit At and k,t{p) and k,l{p) are the variances of 
the transverse and longitudinal momentum transfers per unit At. From the discussion in 
Section UTTbI the precise momentum where fjoip) should be evaluated is known only up to 
corrections of order T/M. 

The Langevin equation is ambiguous until it is discretized. If time is divided into discrete 
steps of At and the the momenta at discrete times are labeled p°, p^, p", then the Ito 
discretization of the Langevin equation is 

where ^'(p") is drawn from the a Gaussian distribution such that 

JTmn 

(e(p")e(P™)> = ^^'(p")^ ' 

and we have defined to coefficients a*(p) and b^^{p) as 

<o(p) = -VDip)p' , (4.6) 
b'^ip) = kl{p) py + kt{p) {6'^ - pY) . (4.7) 

With this choice of discretization, the Langevin equation is equivalent to a Fokker-Planck 
equation, 

The most instructive way to show the equivalence between the Fokker-Plank and Langevin 
approach is to recognize that the Fokker-Planck equation is a Euclidean Schrodinger equation 
and therefore has a phase space path integral representation for the transition probability 
P(p, t|po, to)- (Here the canonical momenta and canonical coordinates are H = ^ and 
Q = p, respectively.) Similarly it is easy to write down a path integral expression for the 
transition probability from the Langevin equations. The two path integrals are the same 
after integrating over canonical momenta |4^. 

Another equally valid discretization is the Stratonovich discretization of the Langevin 
equation 

(pn+i)«_(pn)^ = a;_(p)At + e(p)At. 
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where p = (p"+^ + p")/2. This discretization will lead to a Fokker-Planck equation of the 
a slightly different form 

dP d 1 d dP 

^ + ^(al (p)P) - i-^WJ(p)^) = . 

Clearly the two forms of the Fokker-Planck equation will give the same answer if 

= <1V) - . (4.9) 



The resolution of this ambiguity was clarified by Arnold |42|. The correct procedure is 
to adjust the drag coefficient a* so that the Fokker-Planck equation approaches equilibrium. 
We will discuss the Ito case here. Substituting the equilibrium distribution P oc e~^p^'^ into 
Eq. ()4.8|) and demanding that the r.h.s. return zero we obtain a relation between al^oip) 
and 6*-' (p). This relation, in terms of rjo, kl, and ht is 



(0). ^ _ I^l{v) 



2TE ' 

(4.10) 



Here ci = 3 is the number of dimensions and v is the velocity of the heavy quark. The 
derivative term in square brackets is smaller by a factor of T/ E than the first term and serves 
to renormalize the drag coefficient [i^. The derivative corrections to the drag coefficient 
a* reflect the ambiguity in the momentum at which a* is to be evaluated. Similarly, in the 
Stratonovich case we have r^g'-^* ^ r/g^ + 0{T/E). The 0{T/E) term can be deduced from 
the relation between the Stratonovich and Ito discretizations, Eq. ()4.9p . 

In the previous section we computed in kinetic theory the drag coefficient to leading order 
in T/E, rjj^ . We also computed the longitudinal fluctuations kl to leading order in T/E. 
For the Fokker-Planck approach to be strictly valid the relation i{}^{p) = kl{p)/{2T E) 
must be satisfied; otherwise the stochastic process will not approach equilibrium. This is 
equivalent to the requirement that 

From kinetic theory we have the following equations for the drag and fluctuations: 
f = {fmi±f[k-q'])-f[k-q']{l±fm , 

|<(Ap.)2> = IjMl'ql {f[k]{l±f[k-q'])} . 

In general, these two equations are not related to each other by the fluctuation dissipation 
relation, Eq. ()4.1H) . However, if the transfer energy is small, <^ T, we can expand the 
thermal distribution functions in the drag equation to the leading non-trivial order in u to 
find 

dp 1 



dt 2vT 



\M\'u;' {m{i±fm 

k,q 



|((Ap.)2> = \M\'qnfmi±m)} ■ 
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Upon using the simple kinematic formula u = v Qz given in Appendix El we find the required 
relation Eq. (j4.1H) . Thus only when the transfer energy is small compared to the temperature 
is the Langevin model strictly valid. 

The transfer energy is small only in two limiting cases. In the first case the heavy quark 
is non-relativistic. In this limit the energy transfer u is less than the velocity times the 
momentum transfer, u < vq and therefore is small compared to to the temperature. If the 
quark is relativistic, however, u is only small if the momentum transfer q is small compared 
to T. This is true only to leading-log of T/rriD- Thus we see that the leading-log transport 
coefficients (Eq. (|3.7|) and Eq. (|3.9|) ) satisfy the fluctuation dissipation relation, Eq. (|4.11|) . 

This analysis indicates that an amalgamation of kinetic theory and the Langevin approach 
is correct to leading order. Given the two scales nif) and T, introduce q* between m/j 
and T, say q* ~ y/rriDT. Then treat all collisions with large momentum transfer q > 
q* ^ rriD with ordinary unscreened kinetic theory. The kinetics of the hard collisions 
depends logarithmically on T/q*. All collisions with momentum transfer q < q* T can 
be subsumed into a Langevin process with prescribed transport coefficients that depend 
logarithmically on q*/mo- The dependence on q* cancels when both the Langevin process 
and hard collisions are included. This procedure is only useful when m^) is really much 
smaller than T and we will not adopt it. 

V. A LANGEVIN MODEL FOR HEAVY ION COLLISIONS 

We have seen that, excepting non-relativistic quarks, the Langevin/Fokker-Planck ap- 
proach is a valid description of the kinetics of heavy quarks only to leading logarithm in 
T/tud- Since the coupling is not particularly small, this result says that the full kinetic 
theory should be used to find the evolution of the heavy quark spectrum. 

However, the Langevin model is an appealingly simple framework for studying the 
thermalization of heavy quarks in a heavy ion collision. The Langevin model requires 
two inputs, the drag coefficient ?7^^(f) and the transverse momentum fluctuations kt{v). 
The longitudinal fluctuations are related to ri^^\v) by the fluctuation dissipation relation 
Vo^i"^) ~ i^hiv) / {2TE). (Finally the the drag coefficient can be tweaked as in Eq. ()4.10p 
by terms suppressed by T/E so that the Langevin process approaches equilibrium. ) As 
seen in Fig. |2fb), the longitudinal fluctuations in the Langevin model (^^) are generally 
smaller than the corresponding fluctuations {kl{p)) in the full kinetic theory. The Langevin 
process will underestimate the longitudinal diffusion at high momentum by a factor of 2 to 
4. 

We have employed two models for the drag and transverse momentum diffusion coeffi- 
cients. The first model is based on LO thermal QCD computation of the drag and transverse 
momentum diffusion coefficients. The procedure to set these coefficients is the following. 
First we set rriD/T = 1.5 and then use the results of Section [ill Bl to determine the mo- 
mentum dependence of the coefficients, Vd{p)/vd{0) and kt{p)/ i^t{0). Then ?7d(0) and 
^t(O) = Ki(0) are fixed by specifying the coefficient D and using the relation 



Vd{0) 



1 



T 



2TM 
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This procedure is equivalent to simply adjusting as while leaving mn/T fixed in the equa- 
tions. This model is referred to as LO QCD below. 

The second model is an extreme limit. In this case we set kl equal to a constant, which 
is equivalent to setting 

dp 

— (XV. 

dt 

The transverse and longitudinal diffusion coefficients are then set equal to each other, 
nxiv) = The precise value of kl{v) is specified by setting the diffusion coefficient as 

in the previous model. 

We have also considered a third extreme limit. In this case we set 77d(p) equal to a 
constant, which is equivalent to setting 

dp 

We have found that the results of this model are quite similar to the leading order model 
LO QCD as could have been intuited from Fig. |21^a). We will not discuss this model further. 



A. Solution to the Fokker-Planck Equation for a Bjorken Expansion 



Before considering this Langevin model for the drag and diffusion in detail, let us es- 
timate the effects of thermalization on the spectrum of non-relativistic heavy quarks in 
a medium expanding in a boost invariant fashion. In this section we assume that the 
transport coefficients are evaluated at f = where they are all related by constants, i.e. 
K = kl{0) = k,t{0) = (2TM)?7£i(0). In the non-relativistic limit the Langevin approach is 
strictly valid. 

The Langevin update rule in the local rest frame is 

P 



Ax 
Ap 



M 



:At , 

-r]DP + At 



where ^ is drawn from a Gaussian probability distribution of width (^*^-') 
Boltzmann-Fokker Planck Equation (BFPE) for this update rule is 



dP 



+ 



p 



dP 



dt M dx' 



d_ 

dp^ 



q2 



(5.1) 
(5.2) 

nlAtd'K The 



(5.3) 



where P(x, p, t) denotes the probability to observe a quark with position x and momentum 
p at time t. Our goal in the following paragraphs is to solve this partial differential equation 
for the evolution of the probability distribution for the simple case of a Bjorken expansion. 

For a Bjorken expansion js^, the probability distribution is independent of the trans- 
verse coordinates and invariant under boosts in the z direction. This allows the BFPE to 
be simplified jl^l, and the Green function of the resulting partial differential equation is 
obtained in Appendix ^ To obtain the final spectrum at time t we need to convolve this 
Green function P(p,t|po,to) with the initial spectrum at time 



dN 
d^p dr] 



T]=0,t 



d^PoP{p,t\po,to) 



dN 
d^Po dr] 



(5.4) 



7?=0,io 
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In this restricted sense of Eq. ()5.4p . P(p, t|po, to) is the probabihty that the particle with 
momentum po at time evolves to a particle with momentum p at some time t. P(p, t|po, ^o) 
is given by 



P(p,t|po,to) 



1 



X 



X 



^27c MT^{t) 
1 

1 



exp 



exp 



exp 




2MT^(t) 
2MT^{t) 



(5.5) 



where we have defined the transverse and longitudinal effective temperatures, 

x{t) = [ dt' Mt') , 

J to 

T^{t) = 2 [ dt' T^D{t')T{t')e^^^''^-^^^'\ 



to 



dt' 7]D{t') T{t') 



to 



,2x{t')-2x{t) 



(5.6) 
(5.7) 
(5.8) 



Now we will estimate how the initial spectrum of heavy charm quarks is changed by the 
interactions. For simplicity consider a Bjorken expansion with T(to) = 300 MeV, to = 1 fm, 
M = 1.4 GeV. For an ideal Bjorken expansion, the temperature and drag coefficients follow 



m 

VD{t) 



T{to) 



(5.9) 
(5.10) 



riD^to) is adjusted according to Eq. ()2.8|) to give a specified diffusion coefficient. We will 
compute the spectrum of charm quarks at a final time tf = 6fm. The equilibrium tempera- 
ture at this time is T{tf) = 165 MeV. For the initial spectrum of heavy quarks we will take 
the transverse momentum spectrum from a fit to leading order parton model calculations 
which are described in the next section, 



dN 



drj dyd^pT 



oc S{y-r]) 



(Pt + A 



2\a 



(5.11) 



with a = 3.52 and A = 1.85 GeV. The initial spectrum is shown in Fig. IHland is labeled as 
LO pQCD. Then we convolve this initial condition with the Green function as in Eq. ()5.4|) 
to determine the final spectrum. 

The final transverse momentum spectrum is shown as a function of the diffusion coefficient 
in Fig. ini We observe that the charm spectrum approaches the thermal spectrum only when 
the diffusion coefficient is less than 3/(27rT). This is small relative to the estimates made 
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FIG. 3: (Color online) The transverse momentum spectrum of charm quarks at time tf = 6fm 
for a Bjorken expansion with tq = Ifm and Tq = 300 MeV and T{tf) = 165 MeV. The initial 
transverse momentum spectrum is given by leading order perturbation theory (LO pQCD). 

in the previous section. Further, such a small diffusion coefficient implies a substantial 
suppression of the spectrum at large transverse momentum. This basic observation will be 
quantitatively confirmed when we include radial flow in the next section. As studied in 
Appendix El for large diffusion coefficients interactions simply smear the original spectrum 
(LO pQCD) with a Gaussian. For small diffusion coefficients the spectrum is close to the 
thermal spectrum (T = 165 MeV) up to small viscous corrections. 

B. Elliptic Flow and Suppression of Charm Quarks 

Next we will calculate how the flow of an underlying medium influences the spectrum 
of heavy quarks. If the relaxation time rj'^ is less than the expansion rate of the medium, 
then the heavy quark will follow the medium. If rf'^ is greater than the expansion rate, 
the heavy quark will not follow the medium and the resulting elliptic flow will be small. Of 
course, the relaxation time depends on the momentum of the heavy quark. The goal here is 
to determine the largest possible elliptic flow for a given value of the diffusion coefficient. 

To this end we have placed heavy quarks into a hydrodynamic simulation of the heavy 
ion collision. In the local rest frame of the medium, the heavy quark follows the Langevin 
equations. Further discussion of Lorentz invariance and numerical implementation is given 
below. 

The hydrodynamic simulation is a 2 + 1 boost invariant hydrodynamic model with an 
ideal gas equation of state p = ^e. The temperature is related to the energy density with 
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the Nf = 3 ideal QGP equation of state. We have chosen this extreme equation of state 
because the resuhing radial and elliptic flow are too large relative to data on light hadron 
production. Thus, this equation of state will estimate the largest elliptic flow possible for a 
given diffusion coefficient. 

Aside from the equation of state, the hydrodynamic model is based upon References P. 
At an initial time tq = 1.0 fm, the entropy is distributed in the transverse plane according 
to the distribution of wounded nucleons for a Au-Au collision with an impact parameter of 
b = 6.5 fm. Then one parameter sq = 14fm~^, which is the entropy per unit rapidity per 
wounded nucleon per area, is adjusted to set the initial temperature and total particle yield. 
The value sq = 14fm~^ closely corresponds to the results of full hydrodynamic simulations 
[^, EH and corresponds to a maximum initial temperature of Tq = 265 MeV. 

At the initial Bjorken time tq, the position and momentum distributions of the heavy 
quarks are estimated from leading order parton model calculations. To this end, we have 
distributed the heavy quarks (about a million or so per run) in the transverse plane according 
to the distribution of binary collisions. In the longitudinal direction the heavy quarks are 
distributed uniformly in a large space-time rapidity window t] = —20... 20. Periodic boundary 
conditions are applied in the t] direction. The initial momentum distribution is drawn from 

^^iy-v)7Z^^-TJ^^ ^ (5-12) 



dy drj d'^pr [Px + A^)'^ 

where a = 3.5 and A = 1.849 GeV. This parameterization is a fit to leading order par- 
ton model calculations which were performed with the marvelous CompHEP package [52^ . 
Specifically, charm production at mid-rapidity was computed for proton-proton collisions 
with ^/s = 200 GeV. The charm mass was set to 1.4 GeV and the strong coupling constant 
was evaluated at a scale 4m|. = s [1 — (— t -|- m^/s)^]. Similarly the CETQ4m parton distri- 
butions were evaluated at a scale 4m The resulting momentum distribution is comparable 
to the results of Pythia calculations [5^ . 

A few remarks about Lorentz invariance and numerical implementation are necessary. 
Consider a heavy quark with position and momentum 4-vectors {x')^ and (p')^ in a medium 
with four- velocity u^{x') in the computational frame. Given an infinitesimal time interval 
At', we can calculate Ax' in the computational frame, i.e. Ax' = ^At'. We therefore 
know the 4- vector (Ax')^ and we can compute (Ax)^ in the rest frame of the medium. In 
particular, we know the time interval in the rest frame At. In the rest frame, we can update 
the momentum using the Langevin rule. Now since the quark is on mass shell, the four 
momentum is known in the rest frame and can be boosted back to the computational frame. 
Generally there will be numerical artifacts which are not Lorentz invariant of order At/ L, 
where L is some typical time/length scale of change in the computational frame. Because 
of this dependence on At, it is important to check that Lorentz invariance is respected. 
Several test runs were performed in different frames to verify that numerical artifacts are 
under control at the 1.0% level when the 7 factor is less than 15. The v^At error can 
presumably be avoided by employing a higher order algorithm for stochastic differential 
equations [H^. But, given the complexity of the intermediate boosts, the lowest order 
algorithm was adopted. 

In summary, we place heavy quarks with a reasonable initial transverse momentum dis- 
tribution in a hydrodynamic simulation of heavy ion collisions. Then the heavy quarks are 
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evolved according the the Langevin update rule in the local rest frame. We have simulated 
a solidly mid-central collision, b = 6.5 fm. 

There are at least two items of experimental interest in the heavy ion program. The 
first is the suppression factor Raa, which experimentally is the ratio of the proton-proton 
D meson spectrum to the Au-Au D meson spectrum. We will estimate Raa at mid-rapidity 
using our "input" and "output" charm quark spectrum, 



Raa 



dN 

'^PT J output 



dN 
dpT 



input 



The "input" spectrum has already been described. Given the Langevin nature of the force, 
the spectrum never completely stops changing. However, the drag and fluctuations are 
proportional to powers of the temperature, and therefore their influence decreases at late 
times, as the temperature drops. We have found that the spectrum is nearly frozen after r = 
9 fm. The "output" spectrum is really the spectrum at r = 15 fm. Integrating further would 
have only a negligible effect. A large unknown in the comparison to D meson data is how 
hadronization will affect the final spectrum. In the truly heavy quark limit, hadronization 
would not affect the spectrum significantly. However, the charm quark is not particularly 
heavy, and therefore various models of hadronization will give different answers. Coalescence 
and independent fragmentation are two extreme models which have been studied recently 
[iili^. We will ignore the important effects of hadronization in this work. Therefore, a 



direct comparison of our Raa to the experimental Raa is certainly misguided. The Raa 
factor is illustrated in Fig. HJ^a). This figure is discussed below. 

A second item of experimental interest is elliptic flow. For simplicity we will consider 
only mid-rapidity. Elliptic flow is quantified with V2{pt), 

Id<PdJhz cos(20) 



dydpTdcf) 

V2{Pt) = (cos(20)) = 



y=0 



I' 



dydpxd-<i> 



y=0 



Here the angle is measured with respect to the reaction plane. The experimental deter- 
mination of the reaction plane is discussed in [H^ Again, we will estimate the elliptic 
flow of D mesons with our charm quark spectrum and ignore potential modifications due to 
hadronization. V2{pt) is illustrated in Fig. El^b). 

To illustrate the model dependence, we have also calculated Raa and V2{pt) using two 
extreme models for the drag and fluctuations. As already elaborated, in the first model the 
drag is proportional to the velocity, ^ oc v. The methods and caveats of extracting Raa 
and V2{pt) from the simulation are the same. The results are shown in Fig. El^a) and (b). In 
the second model, the drag is proportional to the momentum, ^ oc p. However, the results 
are not too different from the LO QCD model and therefore the results are not shown. 



C. Comparison with Boltzmann Simulations 

It is useful to compare our Boltzmann-Langevin approach to the classical Boltzmann 
simulations performed by Molnar and subsequent simulations of Zhang et al. ji^. We 
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• D (27tT) = 1 .5 
□ D (27tT) = 3 
A D (27tT) = 6 
O D (27tT) = 12 
T D (27tT) = 24 
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FIG. 4: (Color online) (a) The nuclear modification factor Raa for charm quarks for representa- 
tive values of the diffusion coefficient, (b) V2{pt) for charm quarks for the same set of diffusion 
coefficients given in the legend in (a). In perturbation theory, D x {2'kT) ~ 6 (0.5/as)^. The model 
for the drag and fluctuation coefficients is referred to as LO QCD in the text. The band estimates 
the light hadron elliptic flow for impact parameter h = 6.5 fm using STAR data 




0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 

PT(GeV) PT(GeV) 

FIG. 5: (Color online) (a) The charm quark nuclear modification factor Raa and (b) elliptic flow 
for representative values of the diffusion coefficient given in the legend. In this model, the drag is 
proportional to the velocity, ^ oa v. For further explanation see Fig. ^ 
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will concentrate on Molnar's simulation here. In this simulation the gluon-gluon cross section 
takes the form, ^ oc l/(t — /i^)^, with screening mass, /in = 0.7 GeV. The total gluon-gluon 

cross section fixes as through the relation ctq = The constants are determined by 

comparing the model cross sections at finite t and zero /in with the corresponding unscreened 
perturbative expressions in QCD. The extracted is used to set the charm-quark cross 
sections. To estimate the diffusion coefficient in this model we follow Appendix IB 21 suitably 
modified for the classical statistics in Molnar's simulation. The diffusion coefficient of the 
charm quark in this model can be written in the follow form: 

D= ,^. /fei^i) . (5.13) 



(Cf/Ca) riQCTo \T n. 



where Ug is the density of gluons, Ug is the density of quarks plus anti-quarks, tiq = Ug + Ug 
is the total particle density, and the function f{fio/T,ng/ng) is determined after Appendix 
IB 21 f{fij^/T,ng/ng) is approximately two for typical simulation parameters. 

The magnitude of the resulting elliptic fiow in Molnar's classical Boltzmann simulation, 
with ^ = 2000 for central collisions and ctq = 10 mb, is similar to our results with D ^ 
3/(27rT). Since elliptic fiow develops over a period of l-4fm we will estimate the diffusion 
coefficient at r ^ 2 fm in Molnar's simulation and compare it with D ^ 3/(27rT). The precise 
value of Molnar's diffusion coefficient depends on fi-o/T, which we take to be 2.5. It also 
depends on the ratio of quarks to gluons which we take to be Ug/ug = 36/16 for chemically 
equilibrated 3-fiavor matter and ng/ug = for gluon dominated matter. Using Glauber 
calculations, we estimate the transverse area for central collisions to be A = 100 fm^. Then, 
no = {dN / dy) / [r A) . With these inputs, the diffusion coefficient in Molnar's model is 

(gluon dominated) 

(5.14) 

(chemically equilibrated) 

Taking T 200 MeV and D ^ 3/(27rT), the diffusion coefficient in our model is 

^l..2fm^0.47fm. (5.15) 

Comparing Eq. (j5.14|) and Eq. (j5.15j) we conclude that both approaches agree on the magni- 
tude of the elliptic fiow when the transport coefficients are comparable. Our simulation with 
D = 3/(27rT) corresponds to Molnar's simulation with dN/dy = 2000 for central collisions 
and (To = 10 mb. 




VI. DISCUSSION 

Fig. I^a) and (b) summarize the phenomenological application of our work. Several items 
are apparent in these figures. 

First, we see in Fig. EJ^a) that there is no significant suppression in Raa until the diffusion 
coefficient is of order, D ^ 20/{2ttT). This diffusion coefficient should be compared with 
the perturbative estimate made in Section HTl 

6 /0.5\^ , , 
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The suppression seen in Fig.EJ^a) at large transverse momentum, pt ~ 4.0 GeV, is somewhat 
overestimated since the Langevin description underestimates the longitudinal fluctuations 
by a factor of three to four in this momentum range. However, for px ~ 4.0 GeV radiation 
should become important and further suppress the spectrum. 

Examining Fig. EJ^a) at large momentum, we see that the suppression decreases only 
slightly for px greater than 3 GeV. This is reminiscent of the suppression pattern seen 
in the light hadron data. It seems that all calculation which have a reasonably complete 
model of the fluctuations can reproduce this trend js^ Is^l- This suppression pattern is 
approximately independent of the model for the transport coefficients. Looking at Fig.El^a), 
we see essentially the same pattern when ^ oc f . However, in this case the magnitude of 
the suppression is significantly smaller for a given value of the diffusion coefficient. 

Turning to elliptic flow, we see in Fig. Efb) that the heavy quark elliptic flow rises with 
momentum until ~ 1 — 2 GeV, and then flattens. (When D = 1.5/(27rT), V2{pt) shows a 
qualitatively different behavior which will be discussed shortly.) A similar transition is seen 
in the elliptic flow of light hadrons, illustrated by the band in the figure. For the light hadron 
data, it has been argued that the transition from rising to flat is controlled by transport 
coefficients, and signals a transition from a quasi-thermal to a kinetic regime 0,0- A 
similar argument can be given here. Examining Fig. Ufa), we see that Raa has a funny 
shape at low momentum and then starts to fall exponentially. These features reflect the 
thermal distribution. At higher momentum, Raa stops falling exponentially and approaches 
a constant. This transition coincides with the transition seen in V2{pt) and depends on the 
diffusion coefficient. A detailed study of the correlation between the initial angle of the 
charm quark and the final angle of the charm quark confirms the qualitative conclusion that 
the changes seen in V2 and Raa for ~ 1 — 2 GeV reflect a transition from a quasi-thermal 
to a kinetic regime. Only models which smoothly interpolate between these two regimes can 
effortlessly reproduce these trends in the entire pt range. 

It has been argued that the charm quark could participate in the flow, and thus be 
pushed out to higher transverse momentum ^ i^]. This scenario depends on the diffusion 
coefficient. A detailed analysis of the initial and final momenta of the charm quarks shows 
that, unless the diffusion coefficient is ~ 1.5/(27rT) or less, the charm quarks are not pushed 
to higher momentum by the flow. Indeed, when the diffusion coefficient is greater than 
1.5/(27rT), a particle with final momentum 1.8 GeV typically started with momentum larger 
than this value. When the diffusion coefficient is small, D < 1.5/(27rT), the heavy quarks 
are pushed out to higher momentum by the flow. This push is reflected in the pt dependence 
of the elliptic flow, as illustrated in Fig. Efb). The rapid fall of V2{pt) for pr ~ 3.0 GeV 
indicates that the hydrodynamics is not able to push the heavy quarks beyond this transverse 
momentum for this diffusion coefficient. 

The Pt dependence of the elliptic flow also reflects aspects of the underlying QCD matrix 
elements. In contrast to a constant energy loss model, where ^ oc f , the LO QCD drag 
coefficient rises with the momentum of the heavy quark (see Fig. Efa)). Consequently, in the 
constant energy loss model, the elliptic flow decreases steadily at large momentum rather 
than remain constant as in the LO QCD model. The elliptic flows of the two models are 
compared in Fig. Hfb) and Fig. [Hfb). 

Examining Fig. Efa) and (b), we reach the important conclusion that Raa and 1^2 (pt) 
are tightly correlated. At low momentum, this correlation is encoded in Einstein's relation 



22 



between the energy loss and the hydro dynamic diffusion coefficient, 



dE 
dx 



D 



T 



It is impossible to have a large elliptic ffow without also predicting a significant suppression 
of charm quarks. Indeed, if future measurements find the charm elliptic flow strong, and the 
spectrum not suppressed, the present authors must logically conclude that hydrodynamics 
is not responsible for the observed elliptic flow. This strong conclusion might be mitigated 
if coalescence or some other exotic hadronization mechanism manages to make the D meson 
elliptic flow significantly larger than the underlying c-quark elliptic flow. 

This point is underscored in recent classical Boltzmann simulations by Molnar |4Sj] and 
subsequent simulations by Zhang et al. • We found in Section IV CI that our charm 
quark elliptic flow is comparable to Molnar's results, provided the two simulations have 
approximately the same diffusion coefficient. Our simulation with D = ?>/ (^.-kT) corresponds 
to his simulation with dN/dy = 2000 for central collisions and ctq = 10 mb. However, in 
Molnar's work, the final D-meson elliptic flow is significantly larger than the underlying 
c-quark flow as a result of coalescence. Similarly in the simulation of Zhang et al, the final 
D-meson elliptic flow is amplified by coalescence, although the amplification factor is smaller 
than in Molnar's work due to the details of the coalescence model. Whether coalescence 
can amplify the elliptic flow of D-mesons in a complete dynamical setting remains unclear 



In summary, we have calculated various transport coefficients for a heavy fermion in the 
perturbative Quark Gluon Plasma. Since the gamma factor of the charm quark is not par- 
ticularly large in much of the experimental momentum range, 'jv ^ 4, collisions rather than 
radiation should determine the medium modifications of the heavy quark spectrum. There- 
fore, we have re-examined coUisional energy loss, and also calculated momentum broadening, 
which is essential to determine how the heavy quark spectrum will be modified. To quanti- 
tatively asses the modifications of the heavy quark spectrum, we formulated a Boltzmann- 
Langevin model for the heavy quark kinetics. The model is strictly valid for non-relativistic 
quarks and for relativistic quarks to leading log of T/m^). We characterized the strength 
of the energy loss and momentum diffusion in terms of the diffusion coefficient, which in 
perturbation theory is approximately, D x (27rT) ^ 6(0.5/as)^- The Boltzmann-Langevin 
model was solved analytically for a Bjorken expansion, giving a simple estimate shown in 
Fig. 121 for the medium modifications of the heavy quark spectrum. Then we solved the 
Langevin equations numerically, yielding the modification factor Raa and V2{pt) for charm 
quarks. The results are shown in Fig. |31^a) and (b). 

The diffusion coefficient for a relativistic plasma is of order ~ r^, where is the typical 
collision time of the heavy quark. Thus, future measurements at RHIC will produce a 
tantalizing estimate of this fundamental time scale. Since the diffusion coefficient is a well 
defined concept, this estimate might ultimately be compared to Lattice QCD calculations, 
confirming or rejecting the paradigm of QGP formation in relativistic heavy ion collisions. 
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APPENDIX A: SOLUTION OF THE FOKKER-PLANK EQUATION FORK A 
BJORKEN EXPANSION 



Our starting point is the Boltzmann-Fokker Planck Equation (BFPE) for the probabihty 
P(x, p, t) to find a particle at (x, p), 



dP p' dP 
'dt ^ Mdx^ 



P(x,p,t) . 



(Al) 



It is important to remember that P depends on space and time through the fiow velocity 
and temperature dependences of the medium. 

For a Bjorken expansion, the probability distribution is independent of the transverse 
coordinates and invariant under boosts in the z direction. This allows the l.h.s. of the 
BFPE to be simplified 5l|. Without loss of generality, we can restrict our attention to 
mid-rapidity, z = 0. Demanding that the probability distribution be invariant under small 
longitudinal boosts leads to the requirement that 



dP{z,p^,p, 



dz 



-M 



dP{z,p^,p^ 



2 = 



dpz 



(A2) 



2=0 



where we have used the non-relativistic approximation E = M, and where p-^ denotes the 
transverse vector, i.e. p-^ = {px,Py)- With the standard substitutions ^H, 



Pz 



t 



P{t,pz) = Pit,p, 



and Eq. (jM)), the BFPE reads, 
^ dP 



dP 
'dt 



dpi 



dP 

'dPz 



d^P 



— = SrjnP + r]DPa -^ZI + VdPz— + MTr]D , i + MTr]D [ — 



dpi dpi 



2g2p 



dpi 



(A3) 
(A4) 

(A5) 



This may be turned into the diffusion equation by employing the method of characteristics 
jiij l. With the following substitutions. 



Uz{t) 
Pit.Pi.Pz) 



±_ ff dt ri]j 

Pa 

_ (f dt riD 

p^e-'*o . 



f{t,u^,uz) e'^^o''^- 



(A6) 
(A7) 
(A8) 



Eq. ()A5|) becomes 
dt 



MT e'^^'o + MTr^n ' e^^4 "'^^ . (A9) 



Employing separation of variables and one more change of variables, we write 

f{t,U-^,Uz) = fxiO-^,Ux) fy{9^,Uy) fziOz,Uz) , 



(AlO) 
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with 



= I MTr)D e'^^o"^*"^ 



to 

= I MTrjD (-] 



to 



to 



(All) 
(A12) 



and find three uncoupled diffusion equations, 









dul 


dfy 


_ d'fy 




du^ 




_ 97. 


09, 


dul 



The Green function for the diffusion equation is well known, 

1 



fx{0 ,Ux) 



e 



(A13) 
(A14) 
(A15) 

(A16) 



Unravehng the nested definitions, we find the following Green function for the BFPE for a 
Bjorken expansion: 



Pz{P,t\Po,to) = 



X 



X 



y/2nMT^{t) 
1 

y/27r MT^{t) 

to 1 

t y/27r MT'(t) 



exp 



exp 



exp 



2MT^{t) 

2MT^{t) 

ip" - (t) Po e-^^^^f 
2MT^{t) 



(A17) 



with the intuitive definitions: 



X(t) = [ dt' Tjoit') , 

J to 

T^{t) = 2 f dt' i]D{t')T{t')e^^^''^~''^^'\ 

J to 

T%t) = 2 dt' rjnit') Tit') (j 



2x{t')-2x{t) 



(A18) 
(A19) 
(A20) 



Switching from z to 77 coordinates, we define P(p, i|po,io) = (^/^o)-Pz(p, ^|Po, ^o) and have 



dN 



d^pdrj 



v=o,t 



d^PoP{p,t\po,to) 



dN 



d^Po df] 



(A21) 



ri=0,to 



The differences between the transverse and longitudinal directions of the Green function 
should be noted. 
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Now let us substitute the form expected for the Bjorken solution into these equations. 
Noting that rio cx: T^, we have 



T{t) = T{to) { ^ 



-2a 



(A22) 
(A23) 



Defining xo = Voito) to /{I — 2a) and changing integration variables to integrate over x 
produces 



X = Xo 



l-2a 



1 



T^{t) = T{t) h^{2x + 2xo) - e-'^ T{to) h^{2xo) , 



T^t) = T{t) h^{2x + 2xo) - e-^^ T{to) 

1— 2a 

Here we have defined the function 



/i2^(2xo) 



(A24) 
(A25) 
(A26) 



X 



T e-=^ j^dx' x'^^e""' 7 > -1 
X Ei(a;) 7 = — 1 



and employed the exponential integral Ei(z) = — J^^dt e~^/t. h^{x) may be expressed in 
terms of the incomplete gamma function, but we did not find the results useful in developing 
the series expansions below. For small x we have 



+ ... 



h^{x) 



7 > -1 



(A27) 



7+1 (7+l){7+2) 

(ln(x) + 7b) X - {- \n{x) - 7^; + 1) + ... 7 = -1 
where 7^ fa 0.5772 is the Euler-Mascheroni constant. For large x and all values of 7 we have 

7 , 7(7 - 1) 



h^{x) 



1-- + 

X 



x^ 



+ ... + 0(e-^) . 



(A28) 



For small values of x ^^id xo the heavy quark spectrum is only slightly affected by the 
surrounding medium. Using the small argument expansion of hj{x), we have for small values 
of X and xo , 



T^{t) = T{to) 2r]D{to)to x { 

Similarly, for we have 

T'{t) = T{to) 2rjD{to)to 



to J 



l-3a 



1-3 a 



Mi) 



to 

t I 3 — 3a 



- 1 



a < 



a 



(A29) 



3-3a 



(A30) 
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For small x the effect of rescattering is to convolute the original transverse momentum distri- 
bution with a Gaussian of transverse width MT^{t) and the original longitudinal momentum 
distribution with a Gaussian of width MT^{t). 

For large values of x the system approaches equilibrium and memory of the original 
distribution is lost. The final distribution is a Gaussian which is quite close to the thermal 
distribution. In this limit the effective temperatures in the transverse and longitudinal 
directions are 

Thus, the mean squared transverse momentum is slightly larger than the equilibrium expec- 
tation by an amount of order l/{2ri£,{tf)tf). The longitudinal momentum is slightly smaller 
than the equilibrium value by a similar amount. These corrections are reminiscent of viscous 
corrections to the stress energy tensor. 



APPENDIX B: DETAILS OF LEADING-ORDER CALCULATIONS 

In this appendix we provide the complete details of the derivations of the momentum 
diffusion and energy loss coefficients discussed in the main text. 

1. Matrix elements 

We begin by proving that the matrix elements given in Eq. ()2.9|) are correct. The vacuum 
matrix elements are most easily found in the rest frame of the heavy quark where only t- 
channel gluon exchange contributes. The external gluon polarization sum is performed with 
two transverse polarization vectors which are purely spatial in this frame. The amplitude 
for gluon exchange is ~ P ■ K/ while the amplitude for a Compton process involves 

-Mcon^pton ^ u{pm-^r ~ ^$ - M)-'i'u{p) ~ Tr {-^]f> + ■ (Bl) 

Since is purely temporal, e anti-commutes across it, leading to TWcompton ~ {P ' K)/{P ■ 
K), which is smaller by a factor of k/M than the gluon exchange amplitude. Therefore, the 
Compton amplitude can be ignored. 

Up to a color factor, the squared matrix elements for the scattering of a heavy fermion 
with a quark or gluon via t-channel exchange can be written in the following way: 

1^1' (X L^,iP)M^piK,K')G>'''{Q)G,f3{Q) . (B2) 

Here G^'^{Q) is the internal gluon propagator which we will take to be in Feynman gauge 
ri^'^ IQ^ . Tracing over the heavy fermion line, one easily finds 

L^,{P) = AP^,P, . (B3) 
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This result is independent of the spin of the heavy particle and depends only on the approx- 
imation that its mass is much larger than the of the process. For scattering with a light 
quark, tracing over the light fermion line yields 

4P^P,M^'^ = 16M2A;2(1 + cos d^y) , (B4) 

where k, k' = k — q are the bath particle momenta before and after the collision. For 
scattering with a gauge boson, contracting the three gluon vertex with the polarization 
vectors yields 

A.P^P,M^^ = 4 [P'^iiK^ + K'^)e ■ e' - e,{K - Q) ■ e' - e'^{K' + Q) ■ e)]' . (B5) 
Using e ■ P = e' ■ P = 0, we have 

AP^P^M''" = IQM'^k'^ie ■ ef = 16M^k'^{l + cos^ 4^0 • (B6) 

Next we find a covariant expression for the rest- frame angle cos 6kk'- Since Q is purely 
spatial and K and K' have the same energy in this frame, we have that k' = k — q and 
K''^ = = -2K ■Q + Q2. Further, K ■ K' = k'^{l - cos Okk') = K ■ {K - Q) = -K ■ Q = 
-QV2. Finally, k'^ = {K ■ Pf/M^. Therefore, 

cos^fcfc, = 1-^^^-^. (B7) 

Using this, along with M'^k'^ = {K ■ P)^, allows us to covariantize the expressions we just 
found, resulting in Eq. (j3.6|) . 



Momentum Diffusion — Non-Relativistic Case 



Next complete the calculation of Section ^1 The matrix elements are screened replacing 
the t-channel gluon propagator in Eq. ()B2|) with the HTL propagator. This is most easily 
done in Coulomb gauge, where the HTL propagator in the plasma frame is 



n 



+ 



5, 



ij QiQj 



00 



(B8) 



The substitution —>■ i indicates that only the spatial parts participate. Explicit expressions 
for the HTL self-energies are jH^, IbC^ 



HT(a;,q) 

noo(t^,q) 



m 



m 













+ 










1- 


Yq 



4g3 



In 



In 



q + UJ 
q — uj 



q + uj 
q — UJ 

i 71 



— Z TT 



(B9) 
(BIO) 



Since heavy quark is approximately at rest in the rest frame of the plasma P^ = (M, 0), 
only the temporal part of the gauge boson propagator G^^{Q) enters when computing the 
squared matrix elements in Eq. ()B2|) . For the case of momentum diffusion studied here 
the energy transfer u is small and the HTL correction reduces to simple Debye screening. 
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l/Q^ = l/q"^ 1/(9^ + "^d)- With these observations and the results of the preceding 
section, the screened matrix elements Eq. (j2.9|) are reproduced. 
The total momentum diffusion constant is then, 



3k 



2M 



(27r)3r (k' - q - k)27i6{k' - k)q^ x 



X [Ar^|-M|^.ark%(^)(l-%(A:')) + |A<|J,on^6(fc)(l+n,(A;'))] • (Bll) 

We use the 3-momentum delta function to perform the k' integration and rewrite the the q 
integration as 



(2^ 



47r2 



q dq (icos6'kq 



(B12) 



Using the energy delta function to perform the cos^kq integral and noting the relation, 
cos^kk' = 1 — Q'^/(2/^^), we have 



3k 



47^3 



k^dk 



2k 



qdq 



q- 



Nf 



(g2 + 77^2^)2 

,2 



2 - 



3L 

2fc^ 



X 



'^^k/T _ 1)2 



2 4 

q q 
2 — - I- 

k^ 4fc4 



(B13) 



In performing the q integral, we may drop resulting terms which are subleading in niD/k, 
because the dominant contribution is from k > T, where this ratio is small. The q integral 
in the small m-£)/k limit gives 



3k 



47]-3 



poo 


r 4e 1 




/ k'^dk 


log— -2 




Jo 







N 



ok/T 



f- 



'^k/T + 1)2 



ok/T 



ok/T 



1)2 J 



(B14) 



The remaining k integral is straightforward and yields 



Gtt 



3k 



, 2T 1 

log + - 

mo 2 



lE + 



log + 

rriY) 2 



7. + ^l 



(B15) 

Finally, the relation between k and the diffusion coefficient can be used to determine the 
diffusion coefficient given in the text. 

Finally, let us indicate how Eq. ()B13j) may be used to calculate the diffusion coefficient 
of a heavy quark in Molnar's parton cascade model. In this model, as is related to the 



total gluon-gluon cross section via, (Tq 



47rC|Q| 



where dA = — 1 is the dimension of the 



adjoint representation, and /in is a fixed Debye mass. The particles obey classical statistics 
and therefore /(I ±/) is replaced by Ce^^^'^ , where the constant C is adjusted to reproduce 
the density of the classical particles for each species. With these modifications Eq. ()B13|) 
becomes 



3k 



2T3 



k'^dk 



Ca 



2k 



UqC 



X 



-k/T 



2A;2 



q^ q'^ 



A;2 4A;4 



(B16) 



where Uq is the density of quarks plus anti-quarks and Ug is the density of gluons. The total 
particle density is no = rig + Ug. The form for the diffusion coefficient given in the text 
Eq. ()5.13|) follows from Eq. ()B16|) and the relation between k and D. 
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3. Phase Space 



In order to calculate the transport coefficients we need to integrate over the phase space 
in Eq. ()3.3|) . In this sub appendix we will adapt the integration technology of Baym et al. 



41| to the peculiarities of the heavy quark phase space. 
The desired phase space integration domain is 

We use the spatial delta functions to perform the p' integration and shift variables to in- 
tegrate over the momentum transfer q = k — k' rather than k'. Define the angles cos 9 kq 
and cos 6pq to be the angles between the k and q vectors and the p and q vectors, and 0^;^^ 
to be the angle between the plane containing q,p and the plane containing g, k. Also note 
that p'^ = a/ + (p + q)^- Taking p large, and noting that v = p/p'^ is the velocity, this 
becomes p'^ = p^ + qv cos 6pq plus terms of order q'^/p^ which are small. In the denominator, 
the approximation p'^ = p^ is good enough. 

The phase space becomes, 

r t^hf^ [' d cos 9„ d COS 4, r -p" + k- k') . (B18) 



16(p0)2 (27r)3 7o 
Now we introduce the integration variable uj 



1 = / duj 5{uj - p'^ + p^) , (B19) 



which is the energy transferred to the heavy particle. We now perform the two cos^ inte- 
grations using kinematic relations to rewrite the delta functions as 

5{uj + p-p') = —5 (cos 6'p„ - —) , 
vq \ qv J 

k' ( UJ u? — 

biuj + k' -k) = —si cosOkq - - + H- ^(k - uj) . (B20) 

kq \ q 2kq J 

The phase space integration becomes 



16(p»)M2ir)' i„ V J--±i Jo 2ir 

When q is small, the lower limit on k can be taken to and cosO^q ~ tu/g. 

4. Energy Loss and Momentum Diffusion — Relativistic Case 

Next we will complete the calculation of energy loss and momentum diffusion of Sec- 
tion |^^B] The procedure is straightforward though somewhat cumbersome. We calculate 
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the matrix elements with the HTL propagator, substitute the matrix elements into the ex- 
pressions for the transport coefficients (Eq. (j3.3j) ). and finally integrate over the phase space 
numerically, as parametrized above. 

First we write the matrix elements by inserting the HTL propagator Eq. ()B8|) into 
Eq. ()B2|) and rewrite the resulting expression in terms of the integration variables intro- 
duced in the previous sub appendix. Generally, the matrix elements squared are of the 
form: = A + B cos(0g:fcp) + C cos^ {(pq-kp) ■ After averaging over this azimuthal angle 

the matrix element becomes, {M."^)^ = A + C/2. For a light quark scattering off a heavy 
quark, (A^^)j, becomes 



[cf] I6{p 



,0\2 



Akjk - c^) - [g^ - cu^] [4k{k -uj) + [q^ + u^)] [v^ - - 

2|g2 + nL|^ 4g4|g2_^2 + nT|^ 



(B22) 

where [cf] denotes the color factor, [2CHg'^/2] as in Eq. (j3.6p . 

For gluon scattering, essentially the same procedure is followed. However, the result 
suffers from ambiguities unless the coupling is really small, as discussed in j42|. We will 
follow the prescription of that work and write the gluon matrix element as the scattering 
of a fictitious quark (with the color charge of a gluon) plus an infrared finite remainder. 
To this ephemeral quark we will then apply the HTL corrections and ignore corrections to 
the finite piece. This procedure is correct to leading order since HTL corrections are only 
needed for the eikonal vertex, which is independent of the spin of the hard particle. This 
prescription does something reasonable when m^^/T becomes of order one. Implementing 
this discussion, we compare the quark and gluon matrix elements in Eq. ()3.6|1 and add 

to the quark expression, Eq. ()B22|) . to obtain the gluon squared matrix element, (M^)^ . 
The color factor [cf] is [NcCng"^] for gluon case. We have not written out the M'^/4(P ■ 
term because it is awkward in this choice of integration variables and will be integrated over 
separately. 

To integrate the M^/4(P-i^')^ term we use a different parametrization of the phase space 
integrals. We do not introduce q, but align p with the z axis: 

d^kd^k'd^p' 



= , ^/nN9 / kdkk'dk' [ d cos 9 kpd cos 9 k'p [ - /c' + p° - p °) . 

lQ{p^y {27vy Jo J_i Jq 2-k 

In these variables p'" = p° + vk cos 9kp — vk' cos 9k'p. Next we introduce a new integration 
variable 

/•oo 

1=/ duj6{uj - k{l-v COS 9 kp)) , (B25) 
Jo 

and use the two delta functions to perform the two angular integrations, which yields 



16(p0)2 (27r)3 7o V V Jo 2n 
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The energy transfer is = k' — k and P ■ K = p^k{—l + v cos 6pk) = —ujp^. The vector q is 
comphcated in this coordinate system, but this vector is not needed to integrate M^/4(P • 
K^. With this parametrization, the (pp-kk' integral over M'^/4(P ■ can immediately be 
performed, leaving three integrals that are performed numerically. 

In summary, we take the screened matrix elements insert them into the expressions for 
the transport coefficients and finally perform the phase space integrals. The results of this 
procedure is illustrated in Fig. |21and has been discussed in the main text. 

Analytic expressions for the transport coefficients can be derived to leading logarithm 
in T/m-D [13 ■ We will discuss the energy loss coefficient and leave the details of the other 
transport coefficients to the reader. For q <^ T the matrix elements for the scattering of a 
light quark or gluon on a heavy quark are identical up to a color factor 

Again the color factor [cf] is [2(7/^(7^/2] for quarks and [NcCHg"^] for gluons. For small q and 
uj the weight factor in Eq. (|3.3p is 

^ {n[k-uj]{l ± n[k]) - n[k]{l ± n[k-Lu]) ~ ± ^W) • (B28) 

Next we insert this weight factor and the small q matrix elements into Eq. (j3.3|) and Eq. ()B21|) 
for the energy loss rate. After performing the integral over k, the momentum lass rate is 



dt \ 2 J 12n J 2i;2 V|g2 + noop 2\q'^ - + UtI"^ 

(B29) 

To determine the leading-log coefficient, one should find the coefficient of dq/q when one 
drops the self-energies in this expression. Following this prescription and performing the lj 



integral, we transform Eq. ()B29|1 into 



^ ~ v{N,+ ^] — - — — log log(T/mD) . (B30) 



dt V 2 y 247r 2t;3 1-v 

The coefficient in front of log(T/m£i) is referred to as the leading-log coefficient below. A 
similar analysis but with different weight functions gives the leading log expressions for the 
transverse and longitudinal momentum diffusion coefficients given in the text, Eq. ()3.9p . 

A procedure which is correct to next-to-leading logarithm was formulated by Braaten and 
Yuan [g^I and then used for heavy quark energy loss by Braaten and Thoma [s^. Here the 
momentum integration is divided up into a soft region q < q* and a hard region q* < q, where 
q* is some intermediate scale ttt-d q* '^T . To determine the soft contribution to transport 
coefficient, one should integrate Eq. ()B29|) up to q* with the self energies. For q* ^ mo, 
this integral is a number A^^^iv) plus the leading-log coefficient times log(g*/mD). To 
determine the hard contribution to the transport coefficient, one should drop the self-energies 



in Eq. ()B22jl and numerically integrate the matrix elements over the phase space with q > 
q*. For T ^ q* , this integral is a number A^^^^{v), plus the leading- log coefficient times 
log(T/g*). The sum of the hard and soft contributions is independent of q* and determines 
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V 




Af[v) 




Bf[v) 




Cf{v) 


0.05 


0.0488 


0.7420 


0.0503 


0.7440 


0.0502 


0.7437 


0.10 


0.0567 


0.7499 


0.0617 


0.7570 


0.0620 


0.7568 


0.15 


0.0692 


0.7624 


0.0793 


0.7774 


0.0813 


0.7782 


0.20 


0.0861 


0.7793 


0.1026 


0.8047 


0.1081 


0.8080 


0.25 


0.1074 


0.8006 


0.1312 


0.8383 


0.1425 


0.8464 


0.30 


0.1331 


0.8262 


0.1648 


0.8782 


0.1851 


0.8941 


0.35 


0.1633 


0.8564 


0.2032 


0.9240 


0.2366 


0.9519 


0.40 


0.1981 


0.8913 


0.2462 


0.9758 


0.2981 


1.0211 


0.45 


0.2380 


0.9312 


0.2937 


1.0334 


0.3712 


1.1036 


0.50 


0.2834 


0.9765 


0.3459 


1.0973 


0.4580 


1.2018 


0.55 


0.3348 


1.0280 


0.4029 


1.1677 


0.5617 


1.3195 


0.60 


0.3933 


1.0864 


0.4653 


1.2456 


0.6871 


1.4622 


0.65 


0.4600 


1.1532 


0.5339 


1.3321 


0.8415 


1.6383 


0.70 


0.5371 


1.2303 


0.6101 


1.4293 


1.0369 


1.8616 


0.75 


0.6276 


1.3208 


0.6963 


1.5408 


1.2943 


2.1566 


0.80 


0.7368 


1.4299 


0.7969 


1.6728 


1.6549 


2.5703 


0.85 


0.8747 


1.5678 


0.9200 


1.8374 


2.2119 


3.2097 


0.90 


1.0641 


1.7572 


1.0848 


2.0625 


3.2366 


4.3846 


0.95 


1.3796 


2.0728 


1.3528 


2.4390 


6.0376 


7.5840 



TABLE I: The six functions of velocity A{v), B{v), and C{v) for bosons (6) and fermions (/) 
are defined by Eqs. ()B31MB33|) and parametrize the transport properties of a heavy particle in the 
QGP. 

the transport coefficient. The results of this procedure is the following parametrization of 
the transport properties of a heavy quark in the QGP 



dp 



247r 



log 



2v^ 



1-v 



Nf 



(log(T/mD) + Ab{v)) + ^ (log(T/mD) + Af{v)) 

Cng'T^ /3 1 , {l-v^)\ l+v 

loe 

127r V2 2i;2 Av^ ^ l-v 

N, (log(T/mD) + Bbiv)) + (log(T/mD) + Bf{v)) 



127r 



l-v'^ , l+v 
log 



2v^ 



l-v 



N, (log(r/mD) + C,{v)) + ^ (log(T/mD) + Cj{v)) 



(B31) 



(B32) 



(B33) 



The coefficients A{v), B[v) and C{v) for bosons and fermions are tabulated as a function of 
velocity in Table HI 
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